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TECHNICAL MEMORANDUM 


NASA M ARSHALL ENGINEERING THERMOSPHERE MODEL-VERSION 2.0 

1. BACKGROUND 


The region of the Earth’s atmosphere between about 90- and 500-km altitude is known as 
the thermosphere. The region above -500 km is known as the exosphere. The temperature in the lower 
thermosphere increases rapidly with increasing altitude from a minimum near 90 km but eventually 
becomes altitude independent at upper thermospheric altitudes. This asymptotic temperature, known 
as the exospheric temperature, does not vary with height at any given time due to the extremely short 
thermal conduction time. However, the exospheric temperature does vary with time because of solar 
activity and other factors discussed below. 

The neutral thermosphere is important for two reasons: (1) Even at its low density, it produces 
significant torques and drag on orbiting spacecraft and orbital debris, and (2) the density/height profile 
of the atmosphere above 100- km altitude modulates the llux of trapped radiation encountered at orbital 
altitudes. The state of the neutral thermosphere is most conveniently described in terms of a mean, with 
spatial and temporal variations about that mean. 

The NASA Marshall Engineering Thermosphere (MET) model of the Earth’s atmosphere at 
spacecraft orbital altitudes has evolved, based on work conducted at Marshall Space Flight Center 
(MSEC), over a long pericxl of time. The model is based on the extensive work of Luigi Jacchia and his 
colleagues at the Smithsonian Astrophysical Observatory (SAO) during the 1960’s and 1970’s and was 
developed primarily for application to satellite orbital lifetime prediction, mission analysis relative to 
future reboost requirements, etc. This Technical Memorandum describes the second version of the MET 
model, attributes of the atmosphere at orbital altitudes, and the methods used to model those attributes. 

Short-period atmospheric waves are localized in both time and space in a way that has thus far 
been unpredictable. Consequently, these waves are not explicitly included in NASA Marshall Engineer- 
ing Thermosphere Model- Version 2.0 (MET-V 2.0). Thermospheric winds are also not included. At low 
latitudes (less than 28.5°), wind speeds range from 100 to 200 m/s. At high latitudes (greater than about 
65°), speeds can be 1,500 m/s or more. Rapid (minutes) changes in wind direction (to 180°), probably 
driven by gravity waves, have been observed by a satellite. 



2. DEVELOPMENT 


2.1 Jacchia 1964 Model 

Reference 1 contains a description of the Jacchia 1964 (J64) model. It provides tables of atmo- 
spheric density and composition computed for a wide range of exospheric temperatures, starting from 
a fixed set of boundary conditions at 120 km. The diffusion equation is integrated following empirical 
temperature profiles of exponential form capable of reproducing the densities derived from atmospheric 
drag effects on satellites. Formulae are given that relate the exospheric temperature to solar and geomag- 
netic activity and allow for the diurnal and semiannual variations. The response of the density at the 
200-km level to different types of heating is briefly discussed. 

In the mid-1960’s, personnel of the Aerospace Environment Division at MSEC were responsible 
for providing orbital altitude density inputs for satellite orbital decay prediction. They began studies to 
determine which atmospheric model, when combined with the appropriate orbit propagation program, 
would most accurately predict the observed decay histories of 39 satellites. Since the observed satellite 
decay histories were already available, a number of available models of the thermosphere w'ere selected 
for testing in appropriate orbital lifetime computer programs. Because the satellites had decayed prior to 
the start of the study, actual values of proxy input parameters (described in sec. 3) required by the mod- 
els were used. These proxy parameters were representative of actual solar conditions that occurred 
during the decay periods of the satellites. The J64 model, described in reference 1, had the best perfor- 
mance statistically of the models studied and, therefore, was selected for use by MSEC. 

2.2 Jacchia 1970 Model 

The Jacchia 1970 (J70) model, described in reference 2, is patterned after the 1964 model. The 
main differences consist in the lower height (90 km instead of 1 20 km) of the constant boundary surface 
and a higher ratio of atomic oxygen to molecular oxygen density ([OJ: [O 2 ]) ~ 1.5 at 120 km instead of 
1.0), Mixing is assumed to prevail to a height of 105 km; diffusion occurs above this height. All the 
recognized variations that could be connected with solar, geomagnetic, temporal, and geographic 
parameters are represented by empirical equations. 

Due to the great similarity between the 1964 and 1970 models, and confirmation from observed 
data of the existence of the major differences, the decision was made by MSEC personnel to use the 
newer model. One key difference from the J64 model, in which the temperature-induced density bulge 
remained on the equator all year, was that the J70 model bulge followed the latitudinal excursions of 
the Sun. 
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2.3 Jacchia 1971 Model 


In 1971, SAO published the Jacchia 1971 (J71) model.^ Although an effort was made in the 
earlier models to increase the ratio of [O] to [O 2 ], new observational evidence showed that the increase 
had not been large enough. This model was an attempt to meet as closely as possible the composition 
and density data derived for a height of 150 km by Von Zahn"^ on the basis of all the available mass 
spectrometer and extreme ultraviolet (EUV) absorption data. Mixing is assumed to prevail to a height 
of 100 km; diffusion occurs above this height. All the recognized variations that can be connected with 
solar, geomagnetic, temporal, and geographic parameters are represented by empirical equations. Some 
of the previous Jacchia model equations have been revised, not only in their numerical coefficients but 
also in their form, as a result of new analyses. The J71 model was not considered to be as representative 
of the atmospheric density as the J70 model relative to use in satellite orbital lifetime prediction pro- 
grams. However, it contained some aspects that proved useful in improving MSFC’s thermospheric 
modeling capability. 


2.4 MSFC/J70 Orbital Atmosphere Model 

An unpublished study by environmental personnel at MSFC showed that the J70 model could 
be improved by incorporating the J71 model formulations of seasonal/latitudinal variations in the lower 
thermosphere density below 170 km and seasonal/latitudinal variations in helium (He) above 500 km, 
while retaining the orbital decay prediction accuracy of the J64 model. 

Reference 5 contains a description of the MSFC/J70 Orbital Atmospheric Density model, 
a modified version of the J70 model. The algorithms describing the MSFC/J70 model are provided 
as w'ell as a listing of the computer program. The estimated future 13-mo smoothed values of solar 

flux ( 7 ^ 07 ) ^nd geomagnetic index ( A^), which are required as inputs for the MSFC/J70 model 
when computing future estimates of orbital altitude density conditions, are also discussed. 

2.5 NASA Marshall Engineering Thermosphere Model-1988 Version 

The modifications made and used in the MSFC/J70 Orbital Atmosphere Density model spuri- 
ously introduced step function increases m density requiring minor fairing modifications for their 
elimination. In the late 1980’s, this fairing routine was added, minor programming errors were corrected, 
and the complete program was made more understandable and user friendly. In addition, the calculation 
of atmospheric thermodynamic properties w'as included, and the entire program was updated from 
FORTRAN IV to FORTRAN 77. These modifications and additions produced the first version to be 
known as the MET model used at MSFC. 

Reference 6 contains a description of the MET model-1988 Version, w hich is a modified 
version of the MSFC/J70 Orbital Atmospheric Density model. The modifications to the MSFC/J70 
model required for the MET model are described, graphical and numerical examples of the models are 
included, as is a listing of the MET model computer program. Major differences between the numerical 
output from the MET model and the MSFC/J70 model are discussed. 
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In the course of using this model and plotting the results, it was discovered that density 
discontinuities were occasionally present over certain geographic locations. These were traced to an 
occasional problem with the integration routine used in the model. It w'as replaced with a more accurate 
numerical integration algorithm to solve the problem. Reference 7 contains a detailed description of the 
replacement integration scheme for the integration of the barometric and diffusion equations in the MET 
model. This integration scheme is based upon Gaussian quadrature. Extensive numerical testing reveals 
it to be numerically faster, more accurate, and more reliable than the MSFC/J70 model integration 
scheme (a modified form of Simpson’s rule), which was previously used in the MET model. Numerous 
graphical examples are provided, along with a listing of the modified form of the MET model in which 
subroutine lNTEGRy\TE (using Simpson’s rule) is replaced by subroutine GAUSS (which uses Gaussian 
quadrature). 


2.6 NASA Marshall Engineering Thermosphere Model-1999 Version 

While reviewing the 1988 version of MET, it was detennined that the model would not correctly 
adapt to the change of the approaching new millennium. While performing this update, it was decided 
to also update the calculation of the solar position to use the more recent standard epoch, J2000.0, 
which corresponds to January 1.5, 2000, or Julian date (JD) 2451545.0. This is the standard epoch now 
recommended for use in dynamical astronomy. 

An improvement to the model was also achieved in calculating the solar coordinates. A low 
precision ephemeris for the Sun can be calculated using relatively simple equations.^ This simpler 
scheme for computing the solar position requires fewer arithmetic calculations, yet it maintains more 
than adequate accuracy. The apparent coordinates of the Sun are obtained to a precision of 0°.01; i.e., 
better than an arc minute, and the equation of time to a precision of 0.1 min between 1950 and 2050. 

The other significant modification to the 1988 version of the MET model is in the fairing 
subroutine, which controls the transition from the lower altitude density profile to that of the upper 
altitude region w'here seasonal/latitudinal variation of He is significant. The previous version of this 
subroutine made the transition in a series of small steps. The new scheme introduced in MET-V 2.0 
provides a smooth functional form for the transition. The developmental MET -1999 Version was incor- 
porated into the MSEC Global Reference Atmosphere Model-1999 Version (GRAM-99).^ 
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3. DESCRIPTION 


The MET-V 2.0 is a semiempirical model using the static diffusion method with coefficients 
obtained from satellite drag analyses. It is based on the 1988 version^’^ of MET and work done on the 
1999 version, developed from the Jacchia series of models. Most of the following description also 
applies to the 1988 and 1999 versions of the MET model. With the proper input parameters, specified 
below, an approximate exospheric temperature can be calculated. With e.xospheric temperature specified, 
the temperature can be calculated for any altitude between 90 and 2,500 km from an empirically deter- 
mined temperature profile. In the original development phase of the Jacchia model, the prime objective 
was to model the total neutral mass density of the thermosphere by adjusting temperature profiles until 
agreement between modeled and measured total densities (derived from satellite drag observations) 
was achieved. Agreement between modeled temperatures and those derived from later satellite drag 
and in situ measurements was not always achieved. Thomson-scatter radar temperature measurements 
generally show that the diurnal temperature maximum lags the density maximum by a couple of hours, 
whereas in the MET-V 2.0 model, the temperature and the density maxima and minimal are in phase. 

Studies of the accuracy with which thermospheric models estimate the neutral density have 
shown that an apparent “barrier” exists at the 15-percent level, and models have not thus far been able 
to achieve better performance. A recent study of this type by Marcos et al.^ has shown that, historically, 
this level of accuracy was first achieved by the J70 and J71 models. Figure 1 illustrates this point. The 
MET-V 2.0 model was not explicitly included in this study. However, since MET-V 2.0 incorporates 
the best aspects of Jacchia’s 1970 and 1971 models, it clearly performs at the current 15-percent accu- 
racy state-of-the-art in thermospheric neutral density calculation, while being easier to use than other 
approaches. It was also derived explicitly for use in operational satellite/spacecraft drag applications, 
unlike most thennospheiic density models. The MET-V 2.0 model reflects the J70/J71 models’ accuracy 
results. 
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Figure 1. Statistical evaluation of thermospheric neutral density models. 
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The essence of the MET-V 2.0 model is the calculation of atmospheric density in two major 
regions: the lower thermosphere (altitude 90 km < z < 105 km) and the upper thermosphere (z > 105 km). 
Between the base of the thermosphere (assumed here to be at 90 km) and 105 km, turbulent mixing is 
assumed to predominate, and diffusion dominates at higher altitudes. The density for all points on the 
globe at 90”km altitude is assumed constant and mixing of atmospheric constituents prevails to 105 km. 
Between these two altitudes, the mean molecular mass varies as a result of dissociation of [O 2 ] to [O]. 
An empirical process is employed in the determination of the mean molecular mass distribution between 
90 and 105 km, such that the ratio of [O] to [O 2 ] is 1.5 at 120 km. This makes it agree more closely 
with observations as reviewed by Von Zahn.'^ The input parameters required by the program are altitude, 
latitude, longitude, date (month, day, and year), time (hour and minute), 3 hourly geomagnetic index 
(linear or logarithmic), and the daily 10.7-cm solar radio flux and its average over six solar rotations 
referenced to the midpoint. These parameters and their application will be discussed further in section 4. 
They are essentially the same as provided by the 1988 version of the MET model, but with changes 
made to the computer program to improve the representativeness of the model’s output products 
and efficiency of the computation. 


3.1 Density Distribution 
3.1.1 Barometric and Diffusion Equations 

Density between 90 and 105 km is calculated by integration of the barometric equation. In the 
lower thermosphere, the atmospheric density is computed by integrating the barometric equation: 


Piz) = P90 


Mil) 

f T ) 

''90 

, f 
exp - 

.f 

V y 

h'(z)} 

' V 

•'-V m\z) J 


( 1 ) 


where M is the mean molecular mass, g is the acceleration due to gravity, J’is the local temperature (all 
of w hich are functions of altitude z), and R is the universal gas constant 8.31432 J mole“b The lower 
boundary (ZgQ = 90 km) is assumed to have the following conditions: 

PgQ = 3.46 X 10“^ g cm'"^ 

Ego = 183 K 

Mgg = 28.878 g molerk 

For altitudes above 105 km, the diffusion equation for each of the individual species (O 2 , O, N 2 , 
He, and Ar) is integrated upward from the 105-km level. The number density for an individual species in 
the altitude range 90 < z < 1 05 is calculated using a partition function based upon the sea level composi- 
tion, and this calculation establishes the values for used below. For N 2 , Ar, and He, the number 
density is given by 


«(0 = 


P(z)A^A 

28.96 


( 2 ) 
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where i denotes the species (N,,, Ar, or He), q^(i) is the appropriate volume fraction for species i, and 
is Avogadro’s number. The values of q^(i) used in MET-V 2.0 are 0.781 10 for N 2 , 0.20955 for O-,, 
1.289 X 10”^ for He, and 0.009343 for Ar. For O and O 2 , the equations used are 


n = — =z — — 1 

Miz) [ 


M(z)' 

28.96, 


( 3 ) 


and 


piz)NA\M(z)i, , _ ,1 

n = -= < 1 + qJ 02 ) ■ 

M{z) [28.96^ ^ ^ 


(4) 


respectively. 

In the upper thermosphere (4 > 105 km), the density computation is accomplished by integrating 
the diffusion equation: 


A+a 


A^(Z) = ^105 


LLQ5 

{T(z) 



M{z)g{z) 

RT(Z) 





J 


(5) 


where A is the number density of the species for which the calculation is being done, a is the thermal 
diffusion factor (zero for all species except He, for which a value of -0.380 is used), and the lower 
boundary is at 105 km. It should be noted that, in this case, the molecular mass M is a constant, since 
the concern is for only one species at a time. In both cases, MET-V 2.0 computes the integral using 
Gaussian quadrature. 


The number density of H is assumed to be negligible below' 500 km and in diffusive equilibrium 
above this altitude. Thus, the lower boundary is taken to be at 500 km for this constituent, and its num- 
ber density at this low er boundary is given empirically by 


logloA'500 = 73.13-39.401ogioroo + 5.5(logioToo)^ , 

where is the exospheric temperature given by equation (25). The integration of the diffusion equation 
then proceeds upward from 500-km altitude. 

3.1.2 Composition 

In the lower thermosphere (90 km < z < 105 km), the mean molecular mass is computed 
empirically using the equation 


7 



(7) 


M{ z) '■ 


6 

n=0 


■ 100 )« 


where the coefficients, c^, are used to fit the profile in this region. The result of equation (7) is then 
used in the barometric equation, equation (1), to calculate the total mass density, p. The values of used 
in the MET-V 2.0 model are, following Jacchia,^ 28.15204, -0.085586, 1.284 x lO"'^, -1.0056 x 10“-“’, 
-1.021 X 10-5, ] 0544 X 10-6, ^nd 9.9826 x IQ-*, respectively. 

In the upper thermosphere (z > 105 km), the number densities for the individual species are 
determined by integrating the diffusion equation; i.e., from equation (5), with the values at the lower 
boundary (105 km) given by equations (2)-(4). The total mass density is then given by 


pit)- 


S 

/-I 


M.(7)Ai.(z) 


N 


A 


(8) 


and the mean molecular weight is 


M(z) = 


P(^)A'a 

6 


;=1 


(9) 


where i denotes the individual species (N 2 , Oo, O, Ar, He, or H). 

3.1.3 Gravity Field 

The altitude coordinate system used in the MET-V 2.0 model is referenced to the “surface” 
of the geoid. However, the gravitational field is computed by assuming a spherical planet with 
an effective radius, R-jj = 6356.766 km, where the gravitational acceleration at sea level is assumed 
to be = 9.80665 m s^. These two values are self-consistent and correspond to values of and Rp 
at a latitude of 45° 32' 40" with centrifugal acceleration included. Thus, the gravitational acceleration 
as a function of altitude is given, for purposes of the model, by the equation 

8iz)= 2 ■ ( 10 ) 

1 + -^ 

V 


3.1.4 Temperature Profile 

The temperature varies both temporally and spatially, but its basic profile is defined by a bound- 
ary value TgQ = 183 K and an inflection point at = 125 km. The temperature at the inflection point is 
defined empirically by 
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T^ = a + bT^ + cexp(fcr^) , 


( 11 ) 


where the coefficients a, b, c, and k are 444.3807, 0.02385, 392.8292, and 0.0021357, respectively^ 
and is the exospheric temperature from equation (25). For the low altitude portion of the profile; 
i.e., for 90 < z < 125, the temperature is defined by a fourth degree polynomial: 


4 

T{Z) = T^+Y,<^niz-Zxf , (12 

n=\ 

where the coefficients are determined by requirements placed on the derivatives of the profile. Since 
the lower boundary of the MET-V 2.0 model is taken to be, essentially, the mesopause, the temperature 
there must be T{9Q) = T^q = 183 K, and the gradient there (being a minimum) must be = 0. At the 
inflection point, the second derivative must be zero, and the gradient must satisfy the condition 


4 

— < 
3 


~ ^90 


G,<2 


in order to have no inflections in the region 90< z < z^.. A value of 1.90 for was found 
by experimentation to be best-, so MET-V 2.0 uses 


(13) 
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35 
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Co = 0 


(15) 


C3=-I.7 


W-90V 
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(16) 




(z,.-90)4 ■ 35^ 


(17) 


For z > z^; i.e., z > 125 km, the temperature profile is calculated using the empirical relation 


T(z) = 4- Atan (z - z J 1 4- b(z - z J” | 


(18) 


w 


2 

here zl - — {t^ - tA , B = 4.5 x 10”® for z in km, and n = 2.5. 

K ' 
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3.2 Variations 


3.2.1 Solar Activity 

Solar electromagnetic radiation at ultraviolet (UV) and EUV wavelengths changes substantially 
with the level of solar activity. Thus, thermospheric density is strongly dependent on the level of solar 
activity. An average 1 1-yr solar cycle variation exists; similarly, a “typically” 27-day variation in density 
exists, which is related to the average 27-day solar rotation period. Variations, however, tend to be 
slightly longer than 27 days early in the solar cycle when active regions occur more frequently at higher 
solar latitudes and slightly shorter than 27 days later in the solar cycle when the active regions occur 
more frequently closer to the Sun’s equator. Coronal holes and active longitudes also affect this average 
27-day variation. Changes in the thermospheric density related to changes in level of solar (and geomag- 
netic) activity; e.g., flares, eruptions, coronal mass ejections (CMEs), and coronal holes (CHs), can 
begin almost instantaneously (minutes to hours), although more often a lag of a day or more occurs. 

Various surrogate indices are used to quantify levels of solar activity; one is the 10.7-cm solar 
radio flux, designated /'jo 7 . Although EUV radiation heats the atmosphere, this radiation cannot be 
measured from the ground. The can be measured from the ground and correlates well with the 
EUV radiation. Figure 2 shows the variability of the 7 index during a period of low solar activity. 



Figure 2. Variation in daily F^q 7 flux in 1983. 

The exospheric temperature is computed in the MET-V 2.0 model from a base state temperature 
parameterized using the 10.7-cm solar radio flux. Then several variations due to diurnal variations, 
geomagnetic activity, and semiannual variations are added to the base state to form an empirical estimate 
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of the exospheric temperature. The base stale equation is for the global minimum of the exospheric 
temperature distribution (which occurs at night) when the planetary geomagnetic index (or K^) 
is zero (see sec. 3 for an explanation of these indices): 

T;=383 + 3,32Fio,+1.S(f,„ 7-.F,„7) , d'- 

where T\q 7 is the 10.7-cm (2,800 MHz) solar radio flux in units of Jansky (a Jansky is defined 
as 10“26 W m“2 Hzr^ bandwidth). The bar over Fjq 7 indicates an average over three solar rotations 
prior to the day for which we are doing calculations and three solar rotations afterward; i.e., six solar 
rotations, for a total of 162 days as used by Jacchia in the model development. (See appendix.) 

3.2.2 Diurnal 

Rotation of the Earth induces a diurnal (24-hr period) variation (diurnal tide) in thermospheric 
temperature and density. Due to a lag in response of the thermosphere to the EUV heat source, density 
maximizes around 2 p.m. local solar time for orbital altitudes at a latitude approximately equal to the 
subsolar point. The lag decreases with decreasing altitude. Similarly, minimum density occurs between 
3 and 4 a.m. local solar time at about the negative of the subsolar latitude; i.e., in the diametrically 
opposite hemisphere. In the lowest regions of the thermosphere (120 km and below), where characteris- 
tic thermal conduction time is on the order of a day or more, the diurnal variation is not a predominant 
effect. 


Harmonics of the diurnal tide are also induced in the Earth’s atmosphere. A semidiurnal tide 
(period of 12 hr) and a terdiurnal tide (period of 8 hr) are important in the lower thermosphere (below 
=160 km for the semidiurnal tide and much lower for the terdiurnal tide). Because of large damping 
effects of molecular viscosity, these diurnal harmonic tides are not important at orbital altitudes. 

The diurnal variation is included in MET-V 2.0 under the assumption that the daytime maximum 
Tf^ and nighttime minimum temperatures occur at latitudes which are diametrically opposed to one 
another on the geoid; i.e., they occur at latitudes (p equal to ±5^, respectively, where 5^, is the solar 
declination. This is accomplished by writing the exospheric temperature as 


(l + Rsm^ e) 

1 -1- Rcos” ^ 

cos'” Tj - sin"® 6^ 

\ / 

2 

^ 1 + Rsin'” 9 J 


( 20 ) 


where the fitting parameters are R = 0.31, m = 2.5, and « = 3 and the arguments of the trigonometric 
functions are q = |(|) -i- 5^| / 2, 9 = |(|) -i- 5^| / 2, = H + {5 + p sin (H + Y) . 

The diurnal variation is contained in the expression for in which H is the solar hour angle 
measured from upper culmination (the computation of which will be discussed in section 3.4) and the 
constants j8,p, and /are determined from observations to be -0.6457718, 0.1047198, and 0.7504916, 
respectively; i.e., -37°, 6°, and 43°, respectively. The quantity R is really -1 and is known to be 

correlated with solar activity and with geomagnetic activity.^ Reference 2 provided equations for 
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computing R from annual running means of Kp or from averaged F^q 7 . It was concluded from later 
analyses that the available data were insufficient to establish the rules of the variation, so using a con- 
stant, average value of R was the best approach.^ MET-V 2.0 uses the average value from the J70 
model," as reflected in the 1988 version of MET. 

3.2.3 Geomagnetic Activity 

Interaction of solar wind with the Earth’s magnetosphere (referred to as geomagnetic activity) 
leads to a high-latitude heat and momentum source for the thermospheric gases. Some of this heat and 
momentum is convected to low latitudes. Geomagnetic activity varies, usually having one peak in activ- 
ity just prior to and another just after the peak activity of the solar cycle as defined by the 10.7-cm solar 
radio flux. Also, larger solar cycle peaks are associated with more intense geomagnetic activity. 

A seasonal variation of geomagnetic activity occurs with maxima in March (±1 mo) and September 
(±1 mo) each year. This variation is possibly related to the tilt of the Sun’s rotational axis toward the 
Earth. 


There are two geomagnetic indices that may be used in the MET-V 2.0 model. A logarithmic 
representation of the maximum range over which the magnetic field components vary is designated K. 
Values of K may range from zero to 9, and cover a 3-hr interval. The K values from 12 standard stations 
at latitudes from 48° to 63° and fairly evenly distributed in longitude are averaged to obtain the planetary 
index, K , used in MET-V 2.0. The other index, which may be used in MET-V 2.0, is a , which is a 
linear index derived from and ranges from zero to 400. Sometimes a 24-hr average of the 3 hourly 
values is quoted, in which case the convention is to designate it A^. Although high-latitude ionospheric 
current fluctuations drive the magnetic field fluctuations observed at these stations, the magnetic field 
fluctuations do not drive the thermosphere. Thus, good correlations are not always found between 
observed density changes and the index. Figure 3 shows a year of Op data during a period of low solar 
activity. 



Figure 3. Variation in 3 hourly geomagnetic index, Op, in 1983. 
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When the logarithmic Kp index is used, MET-V 2.0 computes the adjustment to the exospheric 
temperature using 


47’^ = 28 + 0.03 exp(^P , (21) 

but if the linear index is used, the equation 

AT^ = ap+ lOOf 1 - exp(-0.08 rt^)J (22) 

is used. There is a time lag between changes in the geomagnetic indices and temperature variations 
which averages ~6.1 hr; so, whichever index is used, its value should be from a time ~6.7 hr prior 
to the time for which the exospheric temperature is to be computed. 

3.2.4 Semiannual 

This variation in thermospheric density is still poorly understood, but it is believed to be 
a conduction mode of oscillation driven by a semiannual variation in Joule heating in the high-latitude 
thermosphere (as a consequence of a semiannual variation in geomagnetic activity). The variation 
is latitudinal ly independent and is modified by compositional effects. The amplitude of the variation 
is height dependent and variable from year to year with a primary minimum in July, primary maximum 
in October, secondary minimum in January, followed by a secondary maximum in April. Magnitude 
and altitude dependence of the semiannual oscillation vary considerably from one solar cycle 
to another. This variation is important at orbital altitudes. 

The MET-V 2.0 formulation of the semiannual variation is accomplished through the exospheric 
temperature. The variation is roughly proportional to the 10.7-cm solar radio flux and may be computed 
as 


AT = 2.41 + K 


10.7 


0.349 + 0.206sm(2;rT^ + 2260.5)]sm(4;rT^ +247°.6) 


where 


T = - + 0.1145- 
Y 


(2tuI ^ 

-I- sin 1-342®.3 


l2.16 


V r 


J 


with d being the number of days since January' 1 and Y the length of a tropical year; i.e., 365.2422 days. 
Note that d is not the day number of the year. 
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Having computed the various adjustments as described above, the exospheric temperature 
is given by 


'C = T,+A7’^+Ar^ , 

which provides the appropriate value for use in equations (6) and (11) for computing the H number 
density at 500 km and the temperature at the inflection point of the temperature profile, respectively. 

3.2.5 Seasonal/Latitudinal 

The total mass density is modified further by the effects of seasonal/latitudinal density variation 
of the lower thermosphere below 170-km altitude and seasonal/latitudinal variations of He above 
500 km. These two effects were incorporated into MET-V 2.0 using the equations developed by Jacchia 
for his 1971 thermospheric model. ^ 

3.2.5.1 Lower Thermosphere. The barometric and diffusion differential equations were solved 
by assuming fixed boundary conditions at z = 90 km. In reality, the values at that height vary as a func- 
tion of both season and latitude. These seasonal/latitudinal variations are driven in the thermosphere by 
the dynamics of the lower atmosphere (mesosphere and below). Amplitude of the variation maximizes in 
the lower thermosphere between about 105 and 120 km and diminishes to zero around 200 km. Although 
the temperature oscillation amplitude is quite large, corresponding density oscillation amplitude is small. 
This variation is not important at orbital altitudes. 

Computation of this variation in MET-V 2.0 follows the J71 model. ^ In the interest of simplicity, 
the boundary conditions are left fixed, and the density variations are computed using the empirical 
relation 


A\og^Qp(z) = S{z)^Psin 


(26) 


where 


S(z) = 0.014(z - 90)exp -0.0013(z - 90)^ 


( 27 ) 


and 


P 



(j + 100) 


(28) 


This equation is only applied below z = 170 km, since the effect on density variation is negligible above 
this altitude. 
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3.2.S.2 Helium. Satellite mass spectrometers have measured a strong increase in He above the 
winter pole. Over a year, the He number density varies by a factor of 42 at 275 km, 12 at 400 km, and a 
factor of 3 or 4 above 500 km. Formation of the winter He bulge is primarily due to effects of global 
scale winds that blow from the summer to the winter hemisphere. Amplitude of the bulge decreases with 
increasing levels of solar activity due to increased effectiveness of exospheric transport above 500 km, 
carrying He back to the summer hemisphere. Also, a very weak dependence exists of He bulge ampli- 
tude on magnitude of the lower thermospheric eddy diffusivity. 

As was just noted, the seasonal/latitudi nal variation of He decreases with increasing altitude. 
However, He is a sufficiently small fraction of the total density at low altitudes that its seasonal/ 
latitudinal variation is negligible below ~500 km. Above 500 km, MET-V 2.0 uses the equation 


logio «(He) = logio «(He) -I- 0.65 


s. 




sin' 


7t d) 


• 

- sin - 1 
4 


where e is the obliquity of the ecliptie, following the J71 mcxlel.'^ 


(29) 


3.3 Fairing Between Lower and Upper Thermosphere 

In order to prcxluce a smoother, more eontinuous transition from the densities computed below 
500 km to those eomputed above, a fairing algorithm is applied to the base 10 logarithm of the total 
density and He number density for altitudes between 440 and 500 km. The faired density is given by 


p(z) = Piiz)C{z) + P 2 (H [1 - C(z)] , ^30) 

where p^ is the density before adjusting for the seasonal/latitudinal variation of He and p 2 is the density 
afterward, and the fairing coefficient is 

C(z) = cos^ [3(z-Zfc)/2] , (31) 

where z^ is the base fairing height at whieh the fairing caleulation begins, 440 km in this case. The 
fairing coefficients have been chosen such that a gradual transition occurs from the density variation 
below 440 km to the variation above 500 km. 


3.4 Solar Position 

It is evident from the preceding diseussion that a key element in the computation of densities in 
the MET-V 2.0 model is the exospheric temperature, which in turn depends upon the solar hour angle 
and declination. The solar declination is given by 

= sin~^(sinesin/l^) , (3 
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where e is the obliquity of the ecliptic; i.e., the angle between the plane defined by the Earth’s equator 
and the ecliptic plane, and Aj is the solar eeliptic longitude. Due to Earth’s nutation, the obliquity of 
the ecliptic is not constant. (Nutation is the slight pericxlic wobbling motion of the Earth’s rotation axis 
caused by the varying distances and relative directions of the Sun and Moon.) In MET-V 2.0, it is 
computed as 


8 = 23°.439-0°.0000004n , 

where n is the number of Julian days from J2000.0 for the Julian date (JD) corresponding to the date 
specified in the input 


n=JD- 2451 545.0 . (34) 

The Julian date is computed using the Fliegel and Van Flandem algorithm.^® 

The solar ecliptic longitude can be computed by 

=L + l°.915sing + 0°.020sin2g , (^^^ 

where L is the mean longitude corrected for aberrations and g is the mean anomaly; i.e., the mean angle 
from pericenter traveled from some arbitrary starting time assuming constant angular speed. These may 
be computed from 


L = 2800 .460 + 0« .9856474/; 


and 


g = 357° .528 + 0° .9856003/7 , 


(36) 


(37) 


respectively. 

The solar hour angle is the angle between the local meridian and the right ascension of the Sun. 

It is given by 

/y =1-1800.0 + A + /i , (38) 

4 

where t is UT in minutes, A is the longitude of the point for which the calculation is being performed, 
and E is the difference between the apparent solar time and the mean solar time. E is known as the 
equation of time, and it may be computed to a good approximation by 

E = L-a, (39) 

where the solar right ascension a is given by 

a = tan“^(eosetanA,^) . (^*^) 
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3.5 Thermodynamic Quantities 


As an option, MET--V 2.0 computes some thermodynamic quantities, which were computed 
without option in the 1988 version of MET, but were not included in previous neutral thermosphere- 
codes used by MSEC. The change to have these calculations be optional was incorporated to provide 
economy of CPU time for those applieations; e,g., total density for satellite lifetime predictions, that do 
not use these themiodynamic quantities. These include the pressure, pressure scale height, specific heat 
at constant pressure, specific heat at constant volume, and the ratio of the specific heats. The equations 
used to calculate these quantities are presented in sections 3.5.1 through 3.5.3. 

3.5.1 Pressure 

The pressure is calculated directly from the equation of slate under the assumption that the 
atmosphere behaves essentially as an ideal gas. Since the number of moles of gas is given by the ratio 
of the mass of the gas to its molecular weight and the ratio of the mass of the gas to its volume is its 
density, the pressure at altitude z is given by 


pU) = 


Piz) 

M(z) 


RT{7) . 


(41) 


3.5.2 Pressure Scale Height 


The pressure scale height is, by definition, the inverse of the integrand in equation (1). It may be 
observed from equation (1) that this is just the change in height required to change the density of a gas 
of molecular mass M by e~^ in an isothermal atmosphere. To compute this quantity in the MET-V 2.0 
model, it was observed that solving equation (41) for RT and substituting the result into the definition 
of the scale height gives 


H{z) = 


p(z)g(z) 


(42) 


where H(z) is the pressure scale height at altitude z. 

3.5.3 Specific Heat Capacities 

Under the assumption of an ideal gas, it can be showm by statistical mechanical methods that the 
molar specific heat at constant volume, of a monatomic gas is 3i?/2, and for a diatomic gas it is 5R/2. 
It can also be shown from the definitions of and (the molar specific heat at constant pressure) that, 
assuming an ideal gas, is the sum of and the universal gas constant R. Then, the ratio of the specific 
heats, 7 = I , is 5/3 for a monatomic gas and 7/5 for a diatomic gas. Therefore, the ratio of specific 
heats in the thermosphere is calculated in the MET-V 2.0 model using the weighted average 
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( 43 ) 


1.67 n(0) + n(Ar) + n(He) + n(H)] + 1.4 n{0^) + n(N^) 

Y ™ : ± 

n(0) + n(Ar) + n(He) + n(H) + n(02) + n(N 2 ) 
Using the equation of state and the fact that R = Cp - c^ 'm equation (43) yields*^ 


c^Xz) 


M(z)g{z) 

[nz)-i]T(z) 


Then, the specific heat at constant pressure is 


Cp{z)=^y{z)(\Xz) ■ 


(44) 


(45) 
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4. PROGRAM USAGE 


4.1 Input Parameters 

The MET-V 2.0 model (code available from author) is designed to be veiy easy to use. Four 
basic types of input parameters are required by the mtxlel: time, location, solar activity, and geomagnetic 
activity. These may be provided to the model using the sample driver program provided, w'liich reads an 
input data file, puts the data into the INDATA array, and calls the primary subroutine MET. The user is 
free to replace the driver program with his/her own, or to integrate MET-V 2.0 model into a larger 
program, to suit the needs of the project. Note that the INDATA array, and therefore all of its elements, 
is REAL*4. 

Table 1 presents a summary of the contents of the INDATA array. Note that the year is input 
as a four-digit number, while the month, day, hour, and minute are all required to be input as two-digit 
numbers. The first column contains the array element number; the second column identifies the corre- 
sponding parameter; the third column indicates the appropriate range of values, w'here applicable; and 
the fourth column identifies the proper units for that parameter, w-'here applicable. Note that INDATA(9) 
controls whether INDATA(12) is interpreted as or K^, and 1NDATA(13) controls whether the addi- 
tional thermodynamic quantities are calculated. 


Table 1 . INDATA array definition. 


Element No. 

Parameter 

Range 

Units 

1 

Altitude 

90 to 2,500 

Kilometers 

2 

Latitude 

-90 to 90 

Degrees 

3 

Longitude 

-180 to 180 

Degrees 

4 

Year 

1950 to 2050 


5 

Month 

01 to 12 


6 

Day 

01 to 31 


7 

Hour 

00 to 23 


8 

Minute 

00 to 59 


9 

Geomagnetic index type 

1(/(p)or2(ap) 


10 

10.7-cm soiarflux 

0 to 400 

10'* Jansky 

11 

Average 1 0.7-cm soiarflux 

0 to 250 

1 0'* Jansky 

12 

Geomagnetic activity index 

0 to 400 (3p) 
0 to 9 (Kp) 


13 

Thermodynamics flag 

0 ort 
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There are three general categories of applications for the MET-V 2.0 model: (1) After-the-fact 
calculations of densities, (2) real-time (or near real-time) calculations of densities, and (3) future appli- 
cations in which density estimates are required for sometime in the future. In each case, the time and 
location inputs are selected and entered according to the same rules. When entering the time, the year 
is entered in four-digit form, while the month, day, hour, and minute for which the calculation is to be 
done are each entered in two-digit form with the hour and minute being the coordinated universal time 
(UTC). For the position desired, enter the altitude in kilometers and geographic latitude and longitude in 
decimal degrees of the spacecraft location at the time of application. Note that the longitude is measured 
from the Greenwich meridian eastward. The selection of inputs for the solar and geomagnetic activity 
parameters depends upon which application category applies to the calculation as discussed in sections 

4.1.1 through 4.1.3. 

4.1.1 After-the-Fact Calculations 

This type of calculation is required for times more than 81 days in the past, such as the analysis 
of data from a space mission or testing model performance using observed orbital parameters for one 
or more spacecraft. In this case, solar activity is specified using the previous day’s value of the 10.7-cm 
solar radio flux and the centered (about the day for which the calculation is to be done) average of the 
10.7-cm solar radio flux over six solar rotations (162 days), The geomagnetic index, either 
or Kp, is the 3 hourly index value from 6 to 7 hr prior to the time of application. 

4.1.2 Real-Time Calculations 

When the application requires that the density calculations be made for the current time, or less 
than 81 days in the past or future, averaging the observed 10.7-cm flux over six solar rotations as 
described above is clearly not possible. In this case, the previous day’s value of the 10,7-cm solar radio 
flux is still used for the daily F\q 7 value. However, for the average value, observed values are used 
where available, and the estimated future 13-mo smoothed value of the 10.7-cm solar radio flux is 
used as an estimate for the balance of the time period (for which observed daily 7 values are not 
available). These estimated future values of are published by the International Space Environment 
Service (ISES) and are available electronically via the internet (connect to http: //sec.noaa. gov/SolarCycle 
and retrieve the “Table of Predicted Values With Expected Ranges”). The observed a.p or Kp from 
6 to 7 hr prior to time of application is still appropriate for this application category' except for the time 
period up to 81 days in the future, when the estimated f uture value of Ap is used. 

4.1.3 Future Calculations 

When the calculations of densities must be done for a date more than 81 days in the future; e.g., 
for estimating orbital lifetimes, the problem becomes more dependent upon estimates of future average 
values for the 10.7-cm solar radio flux, /'1Q7, and the average planetary geomagnetic index, Ap. Now 
neither the daily nor the 162-day average 10.7-cm solar radio flux is known from observations, nor is the 
3 hourly geomagnetic activity. In this case, the future estimate of the 13 -mo smoothed value of F^q-j 
from ISES, as described in section 4.1.2, is used as an estimate of both the previous day’s value of the 
10.7-cm solar radio flux input and of the 162-day average value. On average, a value of ~12 may be used 
for the 13-mo smoothed estimate of the Ap index for use as an estimate of the 3 hourly geomagnetic 
index input. 
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For future neutral density calculations, the accuracy of the calculation depends primarily upon 
two factors: 

(1) The ability of the thermospheric model to represent the observed neutral density using 
the observed values of solar radio flux (proxy for solar EUV heating) and geomagnetic activity 
used in the development of the model. 

(2) The accuracy with which future solar radio flux and geomagnetic activity can be estimated 
for use as inputs to the thermospheric model. 

As noted in section 3, the MET-V2.0 model has an accuracy of ~15 percent, which is the current 
state-of-the-art for thermospheric density model calculations. However, as noted by Vaughan et al.,^^ 
the major source of uncertainty for future thermospheric neutral density calculations is the estimation 
of future solar EUV heat input. 


4.2 Output Parameters 

Upon completion of the calculation, the MET-V 2.0 model passes the results back to the calling 
routine through the OUTDATA and AUXDATA arrays, both of which are REAL*4. These calculated 
results include the exospheric temperature, temperature at the input altitude (z), number densities for 
each constituent, the mean molecular mass, total mass density and its logarithm, and when desired, 
gravitational acceleration at z, total pressure, pressure scale height, molar specific heat capacities at 
constant pressure and constant volume, and the ratio of the specific heat capacities. All of these param- 
eters are expressed in MKS units. Table 2 details the output arrays. The first column identifies which 
element of which output array contains the data for the parameter in the second column. The notation 
used is AUX for the AUXDATA array (produced when thermcxlynamic quantities are computed) and 
OUT for the OUT DATA array. For example, AUX 2 indicates element 2 of the AUXDATA array, and 
OUT 5 indicates element 5 of the OUTDATA array. The third column gives the MKS units for the 
parameter identified in column 2. The total mass density, temperature, and individual species number 
densities all have the same phase variation in the MET-V 2.0 model. 
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Table 2. Summary of output parameters. 


Location 

Parameter 

Units 

0UT1 

Exospheric temperature 

Kelvin (K) 

OUT 2 

Temperature 

Kelvin (K) 

OUT 3 

N 2 number density 

Molecules per cubic meter (m~3) 

OUT 4 

O 2 number density 

Molecules per cubic meter (m~®) 

OUTS 

0 number density 

Atoms per cubic meter (m"3) 

OUT 6 

Ar number density 

Atoms per cubic meter (m"^) 

OUT 7 

He number density 

Atoms per cubic meter (m"^) 

OUTS 

H number density 

Atoms per cubic meter (m"®) 

OUT 9 

Average molecular weight 

Kilograms per kilo-mole (kg kmoF) 

OUT 10 

Total mass density 

Kilograms per cubic meter (kg 

OUT 11 

Logiotnass density 


OUT 12 

Total pressure 

Pascal (Pa) 

AUX1 

Local gravity acceleration 

Meters per second squared (m s"^) 

AUX2 

Ratio-specific heats 


AUX3 

Pressure scale height 

Meters (m) 

AUX4 

Specific heat constant p 

m2 s“2 K"'i 

AUX5 

Specific heat constant v 

m2 s"2 K"'' 


4.3 Sample Calculation 

Suppose one wants MET-V 2.0 model results, including thermodynamic quantities, for a point 
at lat. 45° N., long. 120° W., and 350-kni altitude on Januaiy 20, 1969, at 19: 1 1 UTC. The 10.7-cm 
solar radio flux for the preceding day (136 x 10^ Jansky) would first be found, then the average 10.7-cm 
solar radio flux be computed over six solar rotations (162 days, ±81 days) centered on January 20, 1969 
(155 X 104.Tansky). Then find the value of the geomagnetic index (rt^ = 9) at «6.7 hr earlier; i.e., 
at about 12:29 UTC on January 20, 1969. Using these values for a sample run with the driver program 
supplied with MET-V 2.0, the input data file (see table 1) reads as follows: 


Altitude 


350.00 


Latitude 


= 45.00 


Longitude 

Year 

Month 

Day 

Hour 

Minute 

Geomagnetic Index Type 
10.7-cm Solar Flux 


= - 120.00 
= 1969.00 
= 1.00 
= 20.00 
= 19.00 

= 11.00 
= 2.00 
= 136.00 


Average 10.7-cm Solar Flux = 155.00 
Geomagnetic Activity Index = 9.00 

Thermodynamics Flag = 1.00 
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Note that in the MET-V 2.0 model longitude is positive toward the east, so 120° west longitude is 
entered as --120. Performing this run wall result in the creation of a file named OlJTPlJT.TXT containing 
the following output (all output in MKS units): 


Exospheric temperature 
Temperature 
N2 number density 
02 number density 
O number density 
A number density 
He number density 
H number density 
Average molecular wt. 
Total mass density 
Log 10 mass density 
Total pressure 
Local grav. acceln. 
Ratio specific heats 
Pressure scale-height 
Specific heat cons, p 
Specific heat cons, v 


= 1031.207 K 
= 1019.849 K 
= 0.3289E+14 /m3 
= 0.1660E+13 /m3 
= 0.281 1E+ 15 /ni3 
= 0.5398E+10 /m3 
= 0.5449E+13 /m3 
= O.lOOOE+07 /m3 
= 17.110 

= 0.9123E-11 kg/m3 
= -11.040 
= 0.452 lE-05 Pa 
= 8.80982 m.sec-2 
= 1.64095 
- 56254.6 m 
= 1244.11 m2, sec-2. K-1 
= 758.168 m2.sec-2.K-l 
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5. CONCLUDING REMARKS 


Using current or past observations of solar radio flux and geomagnetic activity as inputs to 
the MET-V2.0 model will produce thermospheric density estimates with an accuracy of ~15 percent. 
However, using future estimates of these input values from generally accepted statistical models 
(no physical solar model is currently available for use) will result in significantly (order of magnitude 
effects) reduced accuracy for the calculated thermospheric density values. These are key considerations 
to the prediction and statistical confidence of satellite orbital lifetimes, orbital insertion altitudes, reboost 
requirements, etc. for which the MET-V2.0 model (and its predecessors) was developed. 
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APPENDIX --CORRESPONDENCE BY JACK W. SLOWEY, SMITHSONIAN 
ASTROPHYSICAL OBSERVATORY, MARCH 27, 1978 
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SMITHSONIAN INSTITUTION 
ASTROPHYSICAL OBSERVATORY 
60 GARDEN STREET CAMBRIDGE, MASSACHUSETTS 02138 
TELEPHONE 617 864-7910 


March 27, 1978 

D r , W i 1 1 i am W , V a ii g h a n , Chief 
Atmospheric Sciences Division 
ES81, Building 4481 

George C. Marshall Space B’light Center 
Marshall Space Flight Center, Alabama 35812 

Dear Bill: 

It occurs to me that I promdsed to get back to you on the subject of sm^ooth- 
ing to obtain the mean value of the solar flux. This hasn't quite reached the top 
of my list of "things to do today," but I' 11 have a go at it anyway. 

Concerning the gaussian smoothing procedure given on page 20 of the 1977 
Jacchia model, the ” recomonended" value of X should have been 70 days. This is the 
value we've used since we first began to use the method ourselves. Of course, I'm 
sure that we couldn' t distinguish between this and the value of 71 days that is 
quoted. The association with "three solar rotations" may have been the result of 
poor arithmetic following this lesser corruption but, chicken or egg, it is not 
appropriate and should be deleted. 

In the gaussian method, the scale factor X should not be confused with the 
larger interval over which the mean is to be taken. At an interval X from the time 
in question the weight of an included value of F would be 1/e === .37, which is quite 
appreciable. You have to go to an interval of about 3x before the weights really 
become insignificant. It is this limdt, 3x 210 days, that we apply to determine 
"final" values of F with the gaussian n\ethod. It should be emphasized that the 
limit is applied both forward and backward in time: as in the older method, F is 
assumed to be a centered m.ean of the values of F. 

Concerning the older m;ethod, in which F is taken as an unweighted running 
mean, the statement on the sam;e page of the 1977 model that the model 'was based on 
means take-over s total of six solar rotations is correct. The unfortunate as- 
sociation of F with the rnean over three solar rotations that 'was made in both the 
1970 and 1971 models is, again, not a correct reflection of our actual practice. I 
can' t recall the exact chronology, but we may still have been getting F from a 
hand-dra'wn curve through the monthly m.eans when these models were published. In any 
event, when we did finally put the determination of F on a more formal basis, we 
found that a total of six rotations were necessary to represent F "correctly" as 'we 
saiw it. We have used this value (± 3 rotations) ever since. In particular, I have 
used this intervail in connection with this kind of smoothing in all of the model 
computations I' ve ma.de concerning Skylab. 
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Of course, I don' t claim that anything we do with regard to the; de;te;rininaition 
of F is necessarily "best." It is all based on our ow'n subjective feelings in the 
matter. I hope, however, that I have clarified w'hat is consistent with oiir practice 
and "recommended" for use with the atmospheric models that result. 


Sincerely yours. 


i \ 

JaciiJ W. Slawey 


JWS/pl 

?.S. If I can possibly m:eet. deadlines for New York meeting, I' 11 be happy to submit 
a paper. 
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